Advanced Quantitative Methods

Week 6: Panel models

Ozan Aksoy

University of Oxford

Aims for today

  1. Panel models - what and why

    • Research questions
    • Data format

  2. Fixed effects model

    • controlling for time-constant unobserved variables
    • different ways of estimating FE
    • standard errors

  3. Random effects model

  4. Interactions

  5. Dynamic panel models with FE/RE within SEM

  6. Extensions

Research questions

We want to relate changes in X to changes in Y within individuals

For example:

  • how is a change in employment associated with a change in voting behaviour?

  • how are changes in income associated with changes in health?

  • how is a change in parental status associated with changes in employment?

By looking at within individual change, we control for all time-constant variables, observed or unobserved

We usually work with individuals (\(i\)) observed over time (\(t\)), but that is just a specific case of observations nested with some group

Data format

Cases (N) over time (T)

  • small N & mid-large T is usually called a time-series (not covered)
  • large N & T larger than 3 is refered to as panel data

Cases are independent, but observations (T) within a case are correlated

Balanced: all cases have equal T

Unbalanced: T varies from case to case

Regular time intervals

Examples:

  • PSID (US)
  • Understanding Society/BHPS (UK)
  • SOEP (Germany)
  • SHARE (Europe)

Variables varying over \(i\) and \(t\)

time-varying

\(x_{it}\) varies between individuals \(i\) and within individuals across time \(t\)

outcome \(Y_{it}\) also varies between and within cases

time-constant (time-invariant)

\(x_{i.}\) has the same value across \(t\)

this includes fixed characteristics like sex at birth or country of origin

constant in principle vs constant in the data

individual-constant (invariant)

\(x_{.t}\) is the same for all individuals but varies across time

like a level-2 variable in multilevel models, a supra-individual quantity (weather, GDP)

Variance


Some continuous outcome \(y_{it}\)

Grand mean is \(\bar{y}\)

Mean for individual \(i\) is \(\bar{y}_i\)

Total variance \(\sum{(y_{it} - \bar{y})^2}\)

Within variance \(\sum{(y_{it} - \bar{y_i})^2}\)

Between variance \(\sum{(\bar{y_i} - \bar{y})^2}\)


This video has a great visualization of between and within variation in FE models (between minutes 6 and 9)

A simple panel example

id time y x1 x2
1 1 0 0 20
1 2 0 0 21
1 3 6 1 26
2 1 5 1 34
2 2 6 1 36
2 3 7 1 42
3 1 12 1 44
3 2 9 0 30
3 3 9 0 40

N=3, T=3, balanced,
\(y\) is a continuous outcome, \(x1\) binary, \(x2\) continuous

Variance in Y

id time y y_gmean y_imean y_tot y_betw y_with
1 1 0 6 2 -6 -4 -2
1 2 0 6 2 -6 -4 -2
1 3 6 6 2 0 -4 4
2 1 5 6 6 -1 0 -1
2 2 6 6 6 0 0 0
2 3 7 6 6 1 0 1
3 1 12 6 10 6 4 2
3 2 9 6 10 3 4 -1
3 3 9 6 10 3 4 -1

Variance in X1

id time x1 x1_gmean x1_imean x1_betw x1_with
1 1 0 0.556 0.333 -0.222 -0.333
1 2 0 0.556 0.333 -0.222 -0.333
1 3 1 0.556 0.333 -0.222 0.667
2 1 1 0.556 1.000 0.444 0.000
2 2 1 0.556 1.000 0.444 0.000
2 3 1 0.556 1.000 0.444 0.000
3 1 1 0.556 0.333 -0.222 0.667
3 2 0 0.556 0.333 -0.222 -0.333
3 3 0 0.556 0.333 -0.222 -0.333

Variance in X2

id time x2 x2_gmean x2_imean x2_betw x2_with
1 1 20 32.6 22.3 -10.22 -2.33
1 2 21 32.6 22.3 -10.22 -1.33
1 3 26 32.6 22.3 -10.22 3.67
2 1 34 32.6 37.3 4.78 -3.33
2 2 36 32.6 37.3 4.78 -1.33
2 3 42 32.6 37.3 4.78 4.67
3 1 44 32.6 38.0 5.44 6.00
3 2 30 32.6 38.0 5.44 -8.00
3 3 40 32.6 38.0 5.44 2.00

Within model

id time y x1 y_with x1_with x2_with
1 1 0 0 -2 -0.333 -2.33
1 2 0 0 -2 -0.333 -1.33
1 3 6 1 4 0.667 3.67
2 1 5 1 -1 0.000 -3.33
2 2 6 1 0 0.000 -1.33
2 3 7 1 1 0.000 4.67
3 1 12 1 2 0.667 6.00
3 2 9 0 -1 -0.333 -8.00
3 3 9 0 -1 -0.333 2.00

Within model with a time-constant x

id time y x1 y_with x1_with x2_with sexatb sexatb_with
1 1 0 0 -2 -0.333 -2.33 0 0
1 2 0 0 -2 -0.333 -1.33 0 0
1 3 6 1 4 0.667 3.67 0 0
2 1 5 1 -1 0.000 -3.33 1 0
2 2 6 1 0 0.000 -1.33 1 0
2 3 7 1 1 0.000 4.67 1 0
3 1 12 1 2 0.667 6.00 1 0
3 2 9 0 -1 -0.333 -8.00 1 0
3 3 9 0 -1 -0.333 2.00 1 0

No variation in time-constant x
Therefore: can’t explain within-variation in y

Fixed effects panel models in R

  • Two packages: plm and fixest 
## install.packages("plm")
## install.packages("fixest")

  • This video sells the fixest package very well and shows off some of its options

  • plm has handy functions for data inspection and has random effect models

library(plm)
is.pbalanced(df)
[1] TRUE

A within-model by hand

within <- lm(y_with ~ x1_with, data = df)
tidy(within)
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) 2.74e-16     0.282  9.74e-16 1.000   
2 x1_with     4.5 e+ 0     0.732  6.15e+ 0 0.000468

A within-model with “fixest”

library(fixest)
feols(y ~ x1 | id, df)
OLS estimation, Dep. Var.: y
Observations: 9
Fixed-effects: id: 3
Standard-errors: IID 
   Estimate Std. Error t value  Pr(>|t|)    
x1      4.5      0.866     5.2 0.0034782 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.745356     Adj. R2: 0.9375 
                 Within R2: 0.84375

A within-model with “plm”

library(plm)
plm1 <- plm(y ~ x1, data = df, index = c("id"), model="within")
summary(plm1)
Oneway (individual) effect Within Model

Call:
plm(formula = y ~ x1, data = df, model = "within", index = c("id"))

Balanced Panel: n = 3, T = 3, N = 9

Residuals:
   1    2    3    4    5    6    7    8    9 
-0.5 -0.5  1.0 -1.0  0.0  1.0 -1.0  0.5  0.5 

Coefficients:
   Estimate Std. Error t-value Pr(>|t|)   
x1    4.500      0.866     5.2   0.0035 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    32
Residual Sum of Squares: 5
R-Squared:      0.844
Adj. R-Squared: 0.75
F-statistic: 27 on 1 and 5 DF, p-value: 0.00348

A first-differences model

plm2 <- plm(y ~ x1, data = df, index = c("id"), model="fd")
summary(plm2)
Oneway (individual) effect First-Difference Model

Call:
plm(formula = y ~ x1, data = df, model = "fd", index = c("id"))

Balanced Panel: n = 3, T = 3, N = 9
Observations used in estimation: 6

Residuals:
     2      3      5      6      8      9 
-0.833  0.667  0.167  0.167  0.667 -0.833 
attr(,"na.action")
1 4 7 
1 4 7 
attr(,"class")
[1] "omit"

Coefficients:
            Estimate Std. Error t-value Pr(>|t|)   
(Intercept)    0.833      0.312    2.67   0.0557 . 
x1             4.500      0.540    8.33   0.0011 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    42.8
Residual Sum of Squares: 2.33
R-Squared:      0.946
Adj. R-Squared: 0.932
F-statistic: 69.4286 on 1 and 4 DF, p-value: 0.00113

What happens to time-constant variables?

within2 <- lm(y_with ~ x1_with + sexatb_with, data = df)
tidy(within2)
# A tibble: 3 × 5
  term         estimate std.error statistic   p.value
  <chr>           <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  2.74e-16     0.282  9.74e-16  1.000   
2 x1_with      4.5 e+ 0     0.732  6.15e+ 0  0.000468
3 sexatb_with NA           NA     NA        NA       
feols(y ~ x1 + sexatb | id, df)
OLS estimation, Dep. Var.: y
Observations: 9
Fixed-effects: id: 3
Standard-errors: IID 
   Estimate Std. Error t value  Pr(>|t|)    
x1      4.5      0.866     5.2 0.0034782 ** 
... 1 variable was removed because of collinearity (sexatb)
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.745356     Adj. R2: 0.9375 
                 Within R2: 0.84375

Where does our FE come from?

id time y x1 y_with x1_with
1 1 0 0 -2 -0.333
1 2 0 0 -2 -0.333
1 3 6 1 4 0.667
2 1 5 1 -1 0.000
2 2 6 1 0 0.000
2 3 7 1 1 0.000
3 1 12 1 2 0.667
3 2 9 0 -1 -0.333
3 3 9 0 -1 -0.333

Coefficient of 4.5 is identified by cases 1 and 3

One up \(\Delta\) 6 and one down \(\Delta\)-3 make for 9/2=4.5

FE can come from just one case

id time y x1
1 1 0 0
1 2 0 0
1 3 6 1
2 1 5 1
2 2 6 1
2 3 7 1
3 1 12 1
3 2 9 1
3 3 9 1

What do you think the coefficient will be?

FE can come from just one case

feols(y ~ x1 | id, df)
OLS estimation, Dep. Var.: y
Observations: 9
Fixed-effects: id: 3
Standard-errors: IID 
   Estimate Std. Error t value Pr(>|t|)    
x1        6       1.55    3.87 0.011725 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.942809     Adj. R2: 0.9 
                 Within R2: 0.75
  • Fewer cases, but coefficient is bigger
  • Interpretation is now 0 -> 1 (e.g. finding a job / getting marrried)

Changes need not be aligned in time

id time y x1
1 1 0 0
1 2 0 0
1 3 6 1
2 1 5 1
2 2 6 1
2 3 7 1
3 1 12 1
3 2 9 1
3 3 9 0

Changes need not be aligned in time (1)

feols(y ~ x1 | id, df)
OLS estimation, Dep. Var.: y
Observations: 9
Fixed-effects: id: 3
Standard-errors: IID 
   Estimate Std. Error t value Pr(>|t|)    
x1     3.75       1.41    2.66 0.044885 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 1.21335     Adj. R2: 0.834375
                Within R2: 0.585938

Why 3.75?

Case 1: average Y was 0 when x=0 and 6 when x=1, a difference of 6

Case 3: average Y was 9 when x=0 and 10.5 when x=1, a difference of 1.5

Changes need not be aligned in time (2)

id time y x1
1 1 0 0
1 2 0 1
1 3 6 1
2 1 5 1
2 2 6 1
2 3 7 1
3 1 12 1
3 2 9 1
3 3 9 0

Changes need not be aligned in time 2

feols(y ~ x1 | id, df)
OLS estimation, Dep. Var.: y
Observations: 9
Fixed-effects: id: 3
Standard-errors: IID 
   Estimate Std. Error t value Pr(>|t|) 
x1     2.25       1.95    1.16  0.29987 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 1.67498     Adj. R2: 0.684375
                Within R2: 0.210938

Why 2.25?

Case 1: average Y was 3 when x=0 and 6 when x=1, a difference of 3

Case 3: average Y was 9 when x=0 and 10.5 when x=1, a difference of 1.5

Dealing with time: TWFE

Two-way Fixed effects model (TWFE)

Two fixed effects: individual \(i\) and time \(t\)

Our model now looks like this:

\[y_{it} = \beta x_{it} + \alpha_i + \zeta_t + \epsilon_{it}\] time fixed-effects \(\zeta_t\)

indiviudal fixed-effects \(\alpha_i\)

TWFE

In demeaned form TWFE looks like this:

\[(y_{it} - \bar{y_i} - \bar{y_t} + \bar{y}) = \beta (x_{it} - \bar{x_i} - \bar{x_t} + \bar{x}) + (\epsilon_{it} - \epsilon_i - \epsilon_t + \epsilon)\]

TWFE takes into account information about time trends from those without changes in x

TWFE introduces two new assumptions in DiD:

  • parralel trends
  • no anticipation

(see Rüttenauer & Aksoy 2024 for guidance on this)

Are our standard errors ok?

In OLS our errors should be IID: independent and identically distributed

When we have panel data, it is likely that IID doesn’t hold and that our errors are:

  1. autocorrelated: error at t1 is correlated witherror at t2, etc

  2. heteroskedastic: variance is not constant

The normal standard errors, which we would get if we use lm, are not ‘correct’. We are underestimating our se. In panel models it is better to use clustered se that allow for some correlation between errors within individuals.

Works better when you have more clusters, preferable > 50

(for when and which groups to clsuer on, see this paper)

robust standard errors = heteroskadisticy consistent errors, give more weight when residual variance is large

(probably used too often, see Gary King video here)

FE in practice


Enough variance in X?

Effect of going up or down?

  • income: gaining = loosing?
  • marital status: marriage = divorce?

Always:

  • Articulate your theory

  • Look at your data and provide descriptive statistics

  • Think hard about what difference you are identifying and which SE you use

Downsides and caveats of FE

Can not estimate between-differences; no effects of ethnicity, country of origin, etc

If there is very little within variation in a characteristic, that variation can be drowned out by the fixed effects

Coefficients are heavily influenced by cases with a lot of variation on X

Coefficients are weighted averages of ups and downs

Between effects

id time y x1 x2 y_betw x1_betw x2_betw
1 1 0 0 20 -4 -0.222 -10.22
1 2 0 0 21 -4 -0.222 -10.22
1 3 6 1 26 -4 -0.222 -10.22
2 1 5 1 34 0 0.444 4.78
2 2 6 1 36 0 0.444 4.78
2 3 7 1 42 0 0.444 4.78
3 1 12 1 44 4 -0.222 5.44
3 2 9 0 30 4 -0.222 5.44
3 3 9 0 40 4 -0.222 5.44

Between effects

between <- plm(y ~ x1, data = df, index = c("id"), model="between")
summary(between)
Oneway (individual) effect Between Model

Call:
plm(formula = y ~ x1, data = df, model = "between", index = c("id"))

Balanced Panel: n = 3, T = 3, N = 9
Observations used in estimation: 3

Residuals:
        1         2         3 
-4.00e+00  6.66e-16  4.00e+00 

Coefficients:
             Estimate Std. Error t-value Pr(>|t|)
(Intercept)  6.00e+00   6.63e+00     0.9     0.53
x1          -4.26e-15   1.04e+01     0.0     1.00

Total Sum of Squares:    32
Residual Sum of Squares: 32
R-Squared:      0
Adj. R-Squared: -1
F-statistic: 1.68371e-31 on 1 and 1 DF, p-value: 1

Between effects

between <- plm(y ~ x2, data = df, index = c("id"), model="between")
tidy(between)
# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   -7.00      7.03     -0.995   0.502
2 x2             0.399     0.211     1.89    0.309

Between effects

between_ols <- lm(y_betw ~ x2_betw, data = df)
tidy(between_ols)
# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept) 4.71e-16    0.577   8.17e-16 1.000  
2 x2_betw     3.99e- 1    0.0797  5.01e+ 0 0.00155

Remember the within effects?

withinx2 <- plm(y ~ x2, data = df, index = c("id"), model="within")
tidy(withinx2)
# A tibble: 1 × 5
  term  estimate std.error statistic p.value
  <chr>    <dbl>     <dbl>     <dbl>   <dbl>
1 x2       0.301     0.148      2.03  0.0981

Random effects model

mre <- plm(y ~ x2, data = df, index = c("id"), model="random")
tidy(mre)
# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   -4.86      3.97      -1.23 0.220  
2 x2             0.334     0.114      2.94 0.00333

weighted combination of within and between coefficient!

Random effects model

library(lme4)
ml <- lmer(y ~ 1 + x2 +  (1 | id), data = df, REML = FALSE)
tidy(ml)
# A tibble: 4 × 6
  effect   group    term            estimate std.error statistic
  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 fixed    <NA>     (Intercept)       -5.56     3.13       -1.77
2 fixed    <NA>     x2                 0.355    0.0922      3.85
3 ran_pars id       sd__(Intercept)    1.20    NA          NA   
4 ran_pars Residual sd__Observation    1.73    NA          NA

Random effects model in lme4

mwb <- lmer(y ~ 1 + x2_with + x2_imean + (1 | id), data = df, REML = FALSE)
tidy(mwb) # within and between effects, var components unchanged
# A tibble: 5 × 6
  effect   group    term            estimate std.error statistic
  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 fixed    <NA>     (Intercept)       -7.00      4.06      -1.72
2 fixed    <NA>     x2_with            0.301     0.135      2.22
3 fixed    <NA>     x2_imean           0.399     0.122      3.28
4 ran_pars id       sd__(Intercept)    1.16     NA         NA   
5 ran_pars Residual sd__Observation    1.71     NA         NA

Random effects model

RE models takes within and between variance into account

Compare the between and within coefficients in the previous slides

HLM as a panel model

Continuous outcome \(Y\) for individual \(i\) at time points \(j\)

\[Y_{ij} = \beta_{0j} + \beta_{1j}x_{ij} + ... + R_{ij}\]

intercepts can vary by individual:

\[\beta_{0j} = \gamma_{00} + U_{0j}\]

keep coefficients (slopes) constant!

Within individual residuals \(R_{ij}\) have variance \(\sigma^2\)

Between individual residuals have variance: \(var(U_{0j}) = \tau_{00} = \tau_0^2\)

Durbin-Wu-Hausman test

  • Compares FE and RE

  • If estates are close, RE is preferred

  • Better to be led by theory (see here in The Effect)

  • DWH compares FE to basic RE model, other RE models are available

Interactions

  • Can you do interactions in FE models?
    • Yes for time-varying variables, as usual
    • Yes for invariant characteristics
    • Should be led by theory/research question

An interaction example - NLS

ID YEAR gender income nchildren
2.1984 2 1984 2 31500 0
2.1986 2 1986 2 40000 0
2.1988 2 1988 2 20000 0
2.1990 2 1990 2 20100 0
2.1992 2 1992 2 NA 0
2.1994 2 1994 2 NA 1
2.1996 2 1996 2 NA 2
2.1998 2 1998 2 51000 2
2.2000 2 2000 2 45000 2
2.2002 2 2002 2 50450 2

FE model of income and # of children

fe_all <- feols(income ~ nchildren | ID, sample)
summary(fe_all)
OLS estimation, Dep. Var.: income
Observations: 975
Fixed-effects: ID: 95
Standard-errors: IID 
          Estimate Std. Error t value   Pr(>|t|)    
nchildren    25811       3493    7.39 3.4073e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 84,855.3     Adj. R2: 0.288889
                 Within R2: 0.058501

FE by gender

sample$female <- ifelse(sample$gender == 2, 1, 0)
fe_w <- feols(income ~ nchildren | ID, 
              sample %>% filter(female ==1))
tidy(fe_w)
# A tibble: 1 × 5
  term      estimate std.error statistic  p.value
  <chr>        <dbl>     <dbl>     <dbl>    <dbl>
1 nchildren   26738.     3513.      7.61 1.05e-13
fe_m <- feols(income ~ nchildren | ID, 
              sample %>% filter(female ==0))
tidy(fe_m)
# A tibble: 1 × 5
  term      estimate std.error statistic p.value
  <chr>        <dbl>     <dbl>     <dbl>   <dbl>
1 nchildren   23193.     8540.      2.72 0.00703

Quick test

fe_t <- feols(income ~ nchildren + female | ID, sample)
tidy(fe_t)
# A tibble: 1 × 5
  term      estimate std.error statistic  p.value
  <chr>        <dbl>     <dbl>     <dbl>    <dbl>
1 nchildren   25811.     3493.      7.39 3.41e-13

plm comes in handy here to inspect and check the data structure

Interaction

sample$nchildrenXfemale <- sample$nchildren*sample$female
fe_i <- feols(income ~ nchildren + female + nchildrenXfemale | ID, sample)
summary(fe_i)
OLS estimation, Dep. Var.: income
Observations: 975
Fixed-effects: ID: 95
Standard-errors: IID 
                 Estimate Std. Error t value   Pr(>|t|)    
nchildren           23192       6834   3.394 0.00072029 ***
nchildrenXfemale     3545       7952   0.446 0.65580799    
... 1 variable was removed because of collinearity (female)
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 84,845.7     Adj. R2: 0.28824 
                 Within R2: 0.058714

Dynamic panel models with FE/RE within SEM

SEM diagram

Fixed effects as usual

head(d2)
# A tibble: 6 × 5
# Groups:   id [2]
     id      z     t    xit    yit
  <dbl>  <dbl> <int>  <dbl>  <dbl>
1     1 -2.41      1 -3.88  -3.68 
2     1 -2.41      2 -2.13  -3.43 
3     1 -2.41      3 -2.15  -5.13 
4     1 -2.41      4  0.455 -1.96 
5     2  0.555     1 -2.84  -1.25 
6     2  0.555     2 -1.62   0.174
lm_t <- lm(yit ~ xit + z, data = d2)
lm_f <- lm(yit ~ xit , data = d2)
fe   <- plm::plm(yit ~ xit, data = d2, index = c("id", "t"), model = "within")

Summary

modelsummary::msummary(list(lm_t, lm_f, fe), stars = TRUE, 
                       gof_map = c("nobs", "r.squared", "bic") )
(1) (2) (3)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 0.022 -0.027
(0.031) (0.041)
xit 0.416*** 0.636*** 0.418***
(0.014) (0.018) (0.018)
z 0.872***
(0.016)
Num.Obs. 4000 4000 4000
R2 0.560 0.248 0.146
BIC 16816.4 18955.2 15677.8

How to do FE in SEM:

dw <- d2 %>% 
  pivot_wider(names_from = t, values_from = c(xit, yit))

head(dw)
# A tibble: 6 × 10
# Groups:   id [6]
     id      z  xit_1  xit_2 xit_3  xit_4 yit_1   yit_2  yit_3  yit_4
  <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl> <dbl>   <dbl>  <dbl>  <dbl>
1     1 -2.41  -3.88  -2.13  -2.15  0.455 -3.68 -3.43   -5.13  -1.96 
2     2  0.555 -2.84  -1.62  -3.01  1.49  -1.25  0.174  -1.74   1.88 
3     3  2.17  -1.78  -4.57  -2.53 -2.88   5.10 -0.0327  3.49   0.550
4     4 -4.69  -0.177  0.343 -3.32  1.28  -4.96 -7.31   -8.05  -0.255
5     5  0.858  1.89   1.43   1.49  1.83   6.13  0.561   0.729 -1.09 
6     6  1.01   0.344 -2.66  -1.51  3.99   2.73 -0.0258 -0.436 -0.129

Modelling with lavaan 

library(lavaan)
mod <- 'yit_1 ~ b*xit_1 
      yit_2 ~ b*xit_2 
      yit_3 ~ b*xit_3 
      yit_4 ~ b*xit_4 
      A =~ 1*yit_1 + 1*yit_2 + 1*yit_3 + 1*yit_4
      yit_1 +  yit_2 + yit_3 + yit_4 ~ a*1 
      yit_1 ~~ e*yit_1
      yit_2 ~~ e*yit_2
      yit_3 ~~ e*yit_3
      yit_4 ~~ e*yit_4
      xit_1 ~~ NA*A
      xit_2 ~~ NA*A
      xit_3 ~~ NA*A
      xit_4 ~~ NA*A'
fesem <- sem(mod, data = dw)

Results in lavaan 

summary(fesem)
lavaan 0.6-21 ended normally after 41 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        25
  Number of equality constraints                     9

  Number of observations                          1000

Model Test User Model:
                                                      
  Test statistic                               453.934
  Degrees of freedom                                28
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  A =~                                                
    yit_1             1.000                           
    yit_2             1.000                           
    yit_3             1.000                           
    yit_4             1.000                           

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  yit_1 ~                                             
    xit_1      (b)    0.418    0.016   26.929    0.000
  yit_2 ~                                             
    xit_2      (b)    0.418    0.016   26.929    0.000
  yit_3 ~                                             
    xit_3      (b)    0.418    0.016   26.929    0.000
  yit_4 ~                                             
    xit_4      (b)    0.418    0.016   26.929    0.000

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  A ~~                                                
    xit_1             0.817    0.145    5.631    0.000
    xit_2             0.559    0.136    4.098    0.000
    xit_3             0.560    0.138    4.051    0.000
    xit_4             0.574    0.143    4.023    0.000

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .yit_1      (a)   -0.024    0.061   -0.400    0.689
   .yit_2      (a)   -0.024    0.061   -0.400    0.689
   .yit_3      (a)   -0.024    0.061   -0.400    0.689
   .yit_4      (a)   -0.024    0.061   -0.400    0.689
    xit_1             0.090    0.075    1.207    0.228
    xit_2            -0.001    0.072   -0.008    0.994
    xit_3            -0.055    0.073   -0.755    0.450
    xit_4             0.011    0.075    0.144    0.885

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .yit_1      (e)    3.916    0.101   38.730    0.000
   .yit_2      (e)    3.916    0.101   38.730    0.000
   .yit_3      (e)    3.916    0.101   38.730    0.000
   .yit_4      (e)    3.916    0.101   38.730    0.000
    xit_1             5.626    0.252   22.361    0.000
    xit_2             5.149    0.230   22.361    0.000
    xit_3             5.290    0.237   22.361    0.000
    xit_4             5.628    0.252   22.361    0.000
    A                 2.738    0.169   16.211    0.000

Advantages of FE in SEM

  • Error variances can easily be allowed to vary over time
  • FE can have different effects at different times
  • Error auto-correlations can be added
  • Missing data can be handled very easily with FIML (week 8)
  • Other latent variables can be easily added in the model
  • Time-invariant variables can be included in the model (by constraining its correlation with FE to zero)

Disadvantages of FE in SEM

  • Code can be tedious (especially if unbalanced, large T)
  • Convergence issues possible, computationally heavy
  • Does not work well for large T small N

Cross-legged panel model with FE

How to do it in lavaan

mod <- '
 A =~ 1*y2 + 1*y3 + 1*y4
 y2 ~ b*x1 + a*y1
 y3 ~ b*x2 + a*y2
 y4 ~ b*x3 + a*y3
 y2 ~~ x3
 A ~~ y1 + x1 + x2 + x3
 y1 ~~ x1 + x2 + x3
 x1 ~~ x2 + x3
 x2 ~~ x3
 y2 ~~ e*y2
 y3 ~~ e*y3
 y4 ~~ e*y4'
#summary(fit <- sem(mod, data = D))'

Variants of cross-legged models

FE-DPM versus RI-CLPM

FE-DPM

  • Less flexible than Allison et al. (2017)

  • Otherwise very similar

  • Suitable if attention is on x \(\rightarrow\) y

  • But error correlations in Allison et al. makes it more flexible

RI-CLPM

  • Suitable when x and y have stable individual component (psyc construct)

  • Focus is on temporal deviations from stable component

  • Do not block “backdoors” via FE

  • Less “causal” leverage

Further extensions

  1. First differences, between effects, see an example here with interactions (!)

  2. Different types of outcomes: Logit, ordered, multinomial, poisson

  3. Many flavours of standard errors

  4. Combining it with Instrumental Variables

  5. Fixed Effects Individual Slopes (FEIS)

  6. DiD (incl TWFE) models

  7. Present DiD/TWFE as a event study, fect package